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USER'S  MANUAL  FOR  CODES  CHLG  AND  CHTS: 

STEADY  INCOMPRESSIBLE  VISCOUS  FLOW  IN  TWO-DIMENSIONAL  CHANNELS 


^This  report  is  a user's  manual  for  the  two  computer  codes  CHLG  and 
CHTS  for  rapidly  calculating  steady  incompressible  viscous  flow  in  two- 
dimensional  channels^  These  codes  were  both  developed  under  ONR  contract 
1100014-76-C-0324,  administered  under  the  General  Hydrodynamics  Research 
Program  of  the  David  W.  Taylor  Naval  Research  and  Development  Center. 

Copies  of  both  of  these  codes  have  been  sent  to  NSRDC. 

/'Code  CHLG  calculates  laminar  flow  only  in  a channel  with  a general 
shape  for  the  lower  boundary.  Code  CHTS  includes  both  laminar  and  turbulent 
flow  options,  the  latter  by  a simple  eddy  viscosity  model,  but  is  limited  to 
straight  wall  channels.  Both  codes  provide  for  either  Couette  or  Poiseuille 
inflow,  and  use  a symmetry  condition  at  the  top  of  the  mesh. 

This  manual  will  provide  a description  of  input  and  output  variables, 
plus  key  areas  for  possible  modifications  in  the  codes. ^vThe  codes  are  not 
really  user-oriented  software,  but  are  merely  experimen tars codes  used  by 
the  principal  investigator  in  developing  the  semidirect  methods  for  the 
given  problems.  Consequently,  the  codes  will  be  short  on  comment  cards, 
and  will  have  "left-over"  options  from  earlier  experimentations.  Given  the 
limitation  on  time  and  money  under  this  contract,  it  did  not  seem  advisable 
to  expend  a great  deal  of  effort  in  cleaning  up  the  codes. 


Code  CHTS 


Code  CHTS  calculates  laminar  or  turbulent  flow  in  a straight-walled 
channel.  (As  an  option,  it  will  also  calculate  flow  in  a steady  slab  pis- 
ton, but  the  turbulent  flow  option  is  not  operative  for  the  piston  face.) 
The  Input  variables  are  of  two  kinds,  in  DATA  statements  and  in  a READ 
operation.  The  code  Is  designed  to  be  used  on  a time  shared  CDC  6600  com- 
puter system.  At  a certain  point  in  the  program,  the  code  will  print  out 
a request  for  input  variables  which  are  typed  in  at  the  user  terminal.  The 
input  variables  which  appear  in  two  DATA  statements  in  the  main  program 
appear  below.  Those  marked  with  * are  overwritten  by  the  current  form  of 
the  READ  statement  given  later. 
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IL,JL  These  are  the  dimensions  of  the  mesh  in  the  x (stream  wise)  and 
the  y directions.  Currently,  the  variables  are  dimensioned  for 
IL=JL=11.  The  rules  for  changing  dimensions  in  the  main  program 
and  in  Subroutine  EVP9  are  obvious. 

LIM  Maximum  number  of  iterations  allowed,  comnonly  set  to  41. 

NWBC  This  is  the  option  indicator  for  the  vorticity  boundary  condition 

at  the  no-slip  (lower)  wall.  The  suggested  values  are  either 
NVIBC  = 2,  which  gives  Woods'  equation  for  wall  vorticity,  or  NWBC  * 
6,  which  gives  Israeli's  method. 

ELX*,ELY  These  are  the  dimensionless  lengths  of  the  channel  in  the  x and  y 
directions.  Normally,  the  transverse  height  in  y is  the  normaliz- 
ing distance,  so  ELY  = 1.  The  nominal  value  for  ELX  currently  in 
the  program  is  ELX  = 200. 

IPR  This  is  the  option  indicator  for  the  amount  of  printout.  IPR  * 

0 gives  a minimal  printout  including  the  maximum  error  due  to  the 
marching  solution  of  the  linear  equations,  the  converged  values  of 
vorticity  along  the  no-slip  wall,  as  well  as  some  indication  of 
the  wall  vorticity  at  a single  point  and  the  rate  of  change  of  the 
convergence  test  during  each  iteration.  IPR  = 1 gives  a full 
printout  of  the  vorticity  and  stream  function  of  the  converged 
solution,  as  well  as  the  full  errors  in  the  initial  marching  solu- 
tion. IPR  = 2 gives  a great  deal  of  additional  printout,  including 
the  full  influence  coefficient  matrix  used  in  each  of  the  two 
linear  solutions.  It  is  not  recommended  that  IPR  = 2 be  used 
ordinarily. 

ICOR  This  is  the  number  of  iterative  corrections  to  be  used  within  the 
marching  solution  for  the  linear  equation.  It  is  recommended 
that  the  nominal  value  of  ICOR  = 1 be  used. 

NIN  This  option  sets  the  inflow  boundary  conditions  at  I * 1.  NIN  * 1 
gives  Couette  flow,  NIN  * 2 gives  Poiseuille  flow,  for  NBIL  f 1. 
(See  NBIL  below.) 
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NRSTRT  Not  currently  operative. 

NSTACK*  For  NSTACK  » 0,  the  initial  conditions  everywhere  are  set  = inflow 
conditions.  For  NSTACK  * 1,  the  initial  conditions  are  set  from 
the  converged  answer  to  the  previous  case.  For  the  first  case  in 
a RUN,  NSTACK  s 0 is  required. 


CONV* 


This  is  the  convergence  criterion  for  the  non-linear  iterations 
nominally  set  equal  to  5 x 10"4.  For  coarse-mesh  solutions, 

CONV  * 10" 3 Is  entirely  adequate. 

NT  * 0 gives  laminar  flow,  NT  * 1 gives  turbulent  flow  calculated 
with  a simple  eddy  viscosity  model. 

NBIL  t 1 gives  a developing  channel  flow  problem,  while  NBIL  = 1 
gives  a two-dimensional  slab  piston  problem. 

This  is  the  multiplier  for  the  theoretical  under-relaxation  factor 
for  wall  vorticity.  The  recommended  value  is  GFAC  = to  1. 


GILFAC*  This  corresponds  to  GFAC  for  the  piston  face  for  the  option  NBIL  = 
1.  The  recomnended  value  is  GILFAC  = 1. 


NBIL* 


JMDIR 


This  parameter  determines  the  direction  for  the  marching  solution 
for  linear  equation  to  be  in  the  positive  J direction  for  JMDIR  = 
+1.  The  value  JMDIR  = -1  can  also  be  used,  but  the  error  propoga- 
tion  is  better  for  +1. 


When  the  program  Is  compiled  and  run,  the  print  statement  1010  in 
the  code  prints  out  the  label  shown  in  format  statement  905  which  reads: 

"ENTER  RE, ELX, NT, NBIL, NSTACK, GILFAC ;PLEASE" 

The  values  of  the  variables  are  to  be  entered  at  the  user's  terminal.  The 
definitions  follow. 

RE  Reynolds  number  based  on  Inflow  channel  height  and  bulk  mixing 

velocity.  If  RE<0,  the  run  is  terminated  normally. 


ELX  Dimensionless  channel  length  in  the  stream-wise  direction. 

NT  NT  * 0 gives  laminar  flow,  NT  = 2 gives  turbulent  flow. 

NBIL  NBIL  * 0 gives  channel  flow,  and  NBIL  = 1 gives  slab  piston  flow. 

NSTACK  = 0 for  first  case  in  a RUN.  Afterwards,  NSTACK  = 1 sets  initial 

conditions  from  the  converged  answers  of  the  previous  case. 

GILFAC  Multiplier  for  wall  vorticity.  GILFAC  * 1 Is  recommended. 

The  output  of  the  program  comes  in  four  parts.  First  of  all  there 
will  be  an  echo  of  the  input  parameters.  Next,  there  will  be  some  output 
from  the  initial  solution  of  the  linear  equation  by  the  marching  method  In 
Subroutine  EVP9.  The  most  Important  of  these  are  FEMAX  for  IEQ  * 1 and  4, 
which  are  the  maximum  errors  In  stream  function  and  vorticity,  respectively. 
These  should  be  10"4  or  smaller  in  order  to  obtain  reasonable  solutions. 

The  sizes  of  these  quantities  are  strong  functions  of  the  cell  aspect  ratio 
of  the  problem,  which  is  determined  primarily  by  the  parameters  IL,JL,ELX, 
ELY.  A full  description  of  the  error  propagation  characteristics  of  the 
marching  method  will  be  found  in  the  article  "Marching  Methods  for  Elliptic 
Problems,  Part  1",  Volume  1,  No.  1,  of  Numerical  Heat  Transfer  (1978),  pg.  1. 

The  next  set  of  output  is  printed  at  every  iteration  of  the  nonlinear 
method.  The  variables  printed  out  are  the  Iteration  number  NIT,  and  three 
indicators  of  the  convergence  behavior.  The  first  is  DQ,  which  is  the  maxi- 
mum of  the  absolute  value  of  the  change  in  vorticity  since  the  previous 
Iteration.  The  next  value  Is  DQNORN,  which  is  simply  DQ/r  where  r Is  the 
under-relaxation  factor  for  wall  vorticity.  It  Is  this  DQNORM  which  Is 
actually  tested  for  convergence.  It  is  a more  stringent  test  than  testing 
for  DQ,  since  a small  value  of  the  relaxation  factor  could  possibly  give  an 
erroneous  indication  of  convergence.  For  NWBC  * 6,  the  wall  vorticity  is 
not  actually  under-relaxed,  and  DQNORM  = DQ.  The  last  value  printed  out 
at  each  iteration  is  Q(ILC),  which  is  the  actual  value  of  vorticity  along 
the  no-slip  wall  in  the  middle  of  the  mesh  in  the  I direction. 


When  convergence  is  obtained,  the  fourth  set  of  output  is  printed. 

The  values  of  the  mesh  size,  NIT,  CPU  time  in  sec,  and  CPU  time  per  cell 
in  milliseconds  are  printed.  Then,  the  values  of  vorticity  Q along  the 
wall  (that  is,  at  J * 1)  are  always  printed  out  at  the  completion  of  conver- 
gence. If  the  print  option  indicator  IPR  has  been  set  = 1 or  larger,  this 
is  followed  by  a full  field  printout  of  the  variables  stream  function  (S), 
vorticity  (Q),  and  the  velocity  components  U and  V.  The  READ  statement  is 
executed  again,  to  start  another  case;  RE<0  will  terminate  the  run  normally. 

A sample  output  sheet,  including  the  variables  input  by  the  user  from 
the  terminal,  appears  on  page  7. 


Code  CHLG 


Code  CHLG  calculates  laminar  flow  in  a two-dimensional  channel  with 
an  almost  arbitrary  description  of  the  lower  channel  wall.  This  code  re- 
quires the  programs  SGECO  and  SGESL  from  the  UNPACK  Library.  The  input 
variables  of  CHLG  are  very  similar  to  those  for  Code  CHTS,  with  the  following 
exceptions. 

The  variables  VRE,VG,NNRE  and  NNVG  appear  in  DATA  statements.  These 

originally  were  used  to  parameterize  the  Reynolds  number  and  the  under- 

relaxation  factor  at  the  wall  in  DO-Loops.  The  options  are  now  over-written 
by  READ  statements.  The  option  indicator  NMETH  is  set  s 4 in  a DATA  state- 
ment In  order  to  activate  the  Split  NOS  method.  This  is  the  only  option 

that  is  now  active.  Rather  then  set  the  lengths  ELX  and  ELY,  as  in  the 

previous  Code  CHTS,  the  finite  difference  mesh  increments  DX  and  DY  are 
set  in  a DATA  statement,  nominally  DX  * 2 and  DY  * 0.1 . An  additional 
print  option  occurs  In  the  variable  IPRGEO.  For  this  value  * 1 or  greater, 
a detailed  printout  of  the  geometry  of  the  channel  is  given,  including  the 
derivatives  of  the  wall  slope.  The  variable  CELC  which  appears  in  the  first 
DATA  statement  Is  used  to  define  the  length  of  the  channel,  and  will  be  des- 
cribed below.  The  variable  NSOUT  which  appears  in  the  second  DATA  statement 

is  an  outflow  boundary  condition  option.  The  value  NSOUT  a 2 is  recommended. 

/ 

The  READ  statement  requests  the  operator  as  follows: 
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"ENTER  RE , CELC , N JUMP , NWBC , GF , CONV ; PLEASE . (RE<0  TO  STOP)" 


As  before,  RE  is  the  Reynolds  number  based  on  the  inflow  channel 
height,  and  RE<0  terminates  the  RUN  normally.  CELC  defines  the  length 
of  the  channel  according  to  the  relation  ELC  = RE/CELC.  This  makes  the 
length  of  the  channel  to  be  proportional  to  the  Reynolds  number,  allowing 
separated  flow  regions  to  be  maintained  within  the  computational  domain. 

NJUMP  acts  as  NSTACK  in  the  previous  Code  CHTS,  so  NUUMP  * 0 is  required 
for  the  first  case  in  a RUN.  NWBC  is  again  the  boundary  condition  option 
on  wall  vorticity.  Five  different  options  are  still  active  in  the  program, 
but  the  only  recommended  values  are  NWBC  = 2 which  gives  Wood's  equation 
for  vorticity,  NWBC  = 3 which  gives  Jensen's  equation,  and  NWBC  = 6 which 
gives  Israeli's  condition.  The  variable  6F  is  the  multiplier  for  the 
theoretical  under-relaxation  factor  for  wall  vorticity;  the  recommended 
value  is  GF  = 1.  CONV  is  the  convergence  criterion  for  the  nonlinear 
solution;  the  recommended  value  is  CONV  * 10" 3 for  a coarse  mesh. 

The  details  of  the  channel  geometry  are  set  in  Subroutine  GEO.  The 
currently  used  option  is  a channel  that  is  defined  in  Subroutine  CHAN5. 

If  this  option  is  not  used,  the  statement  CALL  CHN5  in  Subroutine  GEO  can 
be  eliminated  in  favor  of  coordinates  for  the  upper  and  lower  surface,  YU 
and  YL  respectively,  being  defined  in  DATA  statements.  Sample  DATA  state- 
ments of  this  type  appear  as  COMMENT  cards  in  Subroutine  GEO.  The  currently 
active  option  in  Subroutine  CHANS  defines  the  lower  wall  of  the  channel 
based  on  a hyperbolic  tangent  function.  The  wall  derivatives  which  are 
necessary  for  the  non-orthogonal  coordinate  transformation  are  evaluated 
numerically  in  Subroutine  GEO.  The  order  of  this  numerical  differentiation 
is  set  by  the  parameter  IODIF,  which  is  set  = 6 at  line  24  of  Subroutine 
GEO.  This  gives  6-th  order  differentiation.  Values  of  IODIF  = 2 or  4 
could  also  be  used. 

With  the  exception  of  the  channel  geometry  data  which  Is  printed 
out  for  IPRGEO  * 1,  the  output  for  Code  CHLG  is  similar  to  the  output  des- 
cribed above  for  Code  CHTS.  Additional  final  output  is  the  value  UWALL  * 
slip  velocity  at  the  lower  wall  (which  is  identically  = 0 in  the  completed 
converged  case)  and  the  % deviation  of  the  vorticity  solution  from  Poiseuille 
flow.  A sample  output  sheet  is  given  on  page  8. 
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FTN  ■>  I -CHTS  » CPT-2 

7.902  CP  SECONDS  COMPILATION  TIME 

✓LGO 

IN  CHTS - IL=  II  JL=  11  L I M=  41  IPR=  0 ICCR-  1 CONV-  .10E-02 
JMDIR*  1 


ENTER : PE  » ELX  • NT  . NB I L « NSTACK « 6 1 LFAC 5 

? 1000. >200. <1 jOjOjI . 


PLEASE 


RE- 10 00.0  ELX  5 ELY-2 0 0 .0  1 . 0 ? DX » DY=  2 0 . 
NWBC-  6 NT=  1 NEIL-  0 UWALL-  0.0  GIL=; 
++*+  TURBULENT  CALCULATION  *+++ 

NWBC-6 « WITH  GW IT  = 2.000 

IN  EVP*  IEC-  4 FEMAX-  .1070E-12 

IN  EVP  > IEC*  1 FEMAX-  .1279E-12 

0 0 .1  0 5 1 

£ 0 . 0 0 0 

=iR* 

£ 0 0 . 0 0 

NIT- 

1 

DG- 

.2000E+0 1 

DCNORM* 

.2000E+01 

o<:  ilc  • l >= 

c! 

. 0 0 0 0 0 0 

NIT- 

u. 

DC- 

. 165 IE* 01 

DGNORM- 

.1651 E+  0 1 

C< ILC  * 1 ) = 

•”i 

.650843 

NIT- 

3 

DG- 

. 1 367E+  0 1 

DCNORM* 

. 1367E+01 

o<:  ilc  > i >= 

5 

.017947 

NIT- 

4 

DG- 

.937 IE* 00 

DCNORM* 

.937 IE* 00 

GX  ILC « 1 >- 

5 

.955027 

NIT- 

5 

DO- 

.4552E+00 

DCNORM* 

.4552E+00 

GK  ILC ? 1 >- 

6 

.410266 

NIT- 

6 

DO- 

. 1656E+00 

DCNORM* 

. 1656E+00 

c<:  ilc  • i > - 

6 

.526933 

NIT- 

% 

DC- 

.5956E-01 

DCNORM* 

•5956E-01 

o<:  ilc  • l > = 

if. 

.520556 

NIT- 

8 

DC- 

.3190E-01 

DCNORM* 

.3190E-01 

c<:  ilc  ? l >= 

I f, 

.488662 

NIT- 

Cl 

DC- 

.3286E-01 

DCNORM* 

•3286E-01 

C< ILC » 1 >- 

6 

.455802 

NIT- 

10 

EG- 

•2613E-01 

DCNORM* 

.2613E-01 

QC ILC 9 1 >= 

if. 

.429673 

NIT- 

11 

DC* 

.1816E-01 

DCNORM* 

.1816E-01 

c<:  ilc  ? i >= 

6 

.411513 

NIT- 

12 

DC- 

. 1 161E-01 

DCNORM* 

. 1 16 IE- 01 

QCILC'l)- 

if. 

. 3999 05 

NIT- 

13 

DC- 

.6957E-02 

DCNORM* 

•6957E-02 

QCILC* 1>- 

if. 

.392948 

NIT- 

14 

DQ- 

.3918E-02 

DCNORM* 

.391SE-02 

CCILC- 1 >= 

if. 

.389030 

NIT- 

15 

DQ- 

.2345E-02 

DCNORM* 

•2345E-02 

C.  I LC?  l :>  = 

6 

.386981 

NIT- 

16 

DC- 

. 1899E-02 

DCNORM* 

. 1899E-02 

C< ILC  » 1 > = 

6 

.386022 

NIT* 

17 

DO* 

.1499E-02 

DCNORM* 

. 1499E-02 

C < ILC  < l ;>  = 

if. 

. 385663 

NIT- 

18 

DC- 

. 1 159E-02 

DCNORM* 

. 1 159E-02 

Q < ILC  » 1 > = 

if. 

.385613 

NIT* 

19 

DQ- 

•8S30E-03 

DCNORM* 

.3S30E-03 

c<:  ILC  * 1 >= 

6 

.385706 

ILjJL 

,3 

11  11 

NIT-  19  CPU  TIME- 

.841  CPU  TIME  PER 

CELL-  8.4 

1 00 


•200000E+01 


.633 089E+01 
.638571E+01 


.633495E+01 

.638571E+01 


.638568E*01 
. 63857 1E+ 01 


. 63857 CE* 01 
. 63857 1E+ 01 


ISEC 

.6 
.6 


38571E+01 
2857 IE* 01 


ENTER:  PE *ELX *NT  jNBIL  j NSTACK  jGILFAC  > PLEASE 

? -1 *1 »1 >1 »1 >1 

.853  CP  SECONDS  EXECUTION  TIME 

/ 
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FTN . I=CHLG  *GPT=2 

11.420  CF  SECONDS  COMPILATION  TIME 

/'LGO 

E NTER  RE  » CELC  » NJUMP  » NWBC  * GF  « CDNV ; PLEASE  . '•  RE<  0 T 0 STOP  > 
'?  1000. » 3 . * 0 ? 6 » 1 . j1  . E-  03 
NIN*  2 NSOUT  = 2 NRST*  40  NSTACK*  0 TCL*  .IE- 03 
CELC*  3.000  ELC=RE "CELC  FOP  SCALING 5 CDNV*  .IE-02 
CHAN?  > XCP.  = .200E+00  F3C  = 10.00 

IN  GEO  t WALL  DERIVATIVES  EVALUATED  BY  I CD  IF*  € 


NPE*  l RE*  . 1 Q00E+04  DX=  .3333E+02  DV=  . 1 Q00E+00  G=  1.0000 

SPLIT  NOS  WALL  EC  OPTION  NO.  6 


NWBC*  6>  WITH  GW IT*  2.000 
PCOND*  .42837 15 IS 12403 IE-01 

IN  EVPIOjIEQ*  4 FEMAX*  .250122E-U 
PCOND*  . 99766896 053368 0E+ 0 0 


IN  EVP 1 0 « I EQ= 

1 FEMAX* 

.234J79E 

-12 

NIT* 

0 

DC* 

.2434E+01 

DONDRM* 

.2434E+01 

0 < ILC  « 1 > = 

.731952 

NIT* 

1 

DO* 

. 1 1 08E+01 

DQNDRM* 

. 1 1 08E+ 0 1 

0 < I LC  » 1 > = 

-.107800 

NIT* 

2 

DO* 

•5451E+00 

DONGRM* 

.545 1E+00 

0< ILC  * 1 > = 

-.177841 

NIT* 

•Z( 

DO* 

.4359E+00 

DQNDRM* 

.4559E+00 

Q< ILC  ? 1 >* 

. 020898 

NIT* 

4 

DO* 

. 154£E*00 

DONGRM* 

. 1542E+00 

0 < ILC  • 1 > = 

. 085050 

NIT* 

5 

DO* 

. 1 059E+00 

DQNDRM* 

. 1 059E+  0 0 

o<:  ilc  » i >= 

. 154033 

NIT* 

6 

DO* 

.527 IE- 01 

DONGRM* 

.527 IE- 01 

Q<ILC»1>= 

.149927 

NIT* 

DO* 

.3234E-01 

DONGRM* 

.3234E-01 

0<ILC*1>* 

.131724 

NIT* 

3 

DO* 

. 1S8SE-01 

DONORM* 

. 1885E-01 

o<:  ilc  . i >* 

.126621 

NIT* 

9 

DO* 

.8448E-02 

DQNDRM* 

.S44SE-02 

Qc:iLC«l>  = 

• 1 2825 1 

NIT* 

10 

DO* 

•5298E-02 

DONORM* 

.5298E-02 

Q<ILCf 1)= 

.130656 

NIT* 

11 

DO* 

•3376E— 02 

DONORM* 

•3376E-02 

o<:  ilc  > i >* 

.131708 

NIT* 

12 

DO* 

. 1423E-02 

DONORM* 

. 1423E-02 

Q < 1 LC « 1 > = 

. 131479 

NIT* 

13 

DO* 

.9121E-03 

DQNDRM* 

.9 12 IE- 03 

Q<ILC»1>= 

.130962 

IL  > JL=  11  11  NIT*  13  CPU  TIME  = .724  CPU  TIME  PER  CELL  = 7.240  MSEC 

QWALL  AND  UWALL  FOLLOW 

.300000E+01  . 152359E+01  .143376E+00  -.115882E+00  -.442829E-01  .130962E+00 
.246201E+00  .382672E+00  .425282E+00  .533015E+00  .537695E-r00 


0. 


. 976383E- 05  - . 28766 0E- 03  - . 1 247? 0E-03 
-.933845E-04  -.647472E-0 4 -.686358E-04 


-.14  0475E- 03  - . 1 2 0224E- 03 
-.490483E-05  0. 


7.  DEVIATION  FROM  FULLY  DEVELOPED  PGISEUILLE  FLOW 
0.  - .374 026E+02  -.392318E+02  - . 1 13795E+03 

- .659929E+02  - .466594E+02  - .4 019? IE* 02 


- . 1 059  0*E+ 03  - .82 1 07SE+  02 
- . 24387 1E+ 02  - .243871E*02 


ENTER  PE •CELT  .NJUMP -NWBC jGF jCONV ! PLEASE . < RE< 0 TO  STOP > 

? -1  »1  «1  *1  .1  >1 

.7 48  CF  SECONDS  EXECUTION  TIME 


